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ABSTRACT 

The  transient  analysis  of  hydrodyr,am' c  lubrication  of  a  point-contact  is 
presented.  A  body-fitted  coordinate  system  is  introduced  to  transform  the 
physical  domain  to  a  rectangular  computational  domain,  enabling  the  use  of  the 
Newton-Raphson  method  for  determining  pressures  and  locating  the  cavitation 
boundary,  where  the  Reynolds  boundary  condition  is  specified.  In  order  to 
obtain  the  transient  solution,  an  explicit  Euler  method  is  used  to  effect  a 
time  march.  The  transient  dynamic  load  is  a  sinusoidal  function  of  time  with 
frequency,  fractional  loading,  and  mean  load  as  parameters. 

Results  Include  the  variation  of  the  minimum  film  thickness  and  phase-lag 
with  time  as  functions  of  excitation  frequency.  The  results  are  compared  with 
the  analytic  solution  to  the  transient  step  bearing  problem  with  the  same 
dynamic  loading  function.  The  similarities  of  the  results  suggest  an 
approximate  model  of  the  point  contact  minimum  film  thickness  solution. 


*NASA  Resident  Research  Associate  at  Lewis  Research  Center. 


INTRODUCTION 


Nonconformal  contact  machine  elements  in  power  train  systems  such  as 
gears,  rolling  element  bearings,  and  cam  and  follower  mechanisms  are  subject 
to  transient  lubrication.  The  transient  characteristics  are  due  to  the  time 
variation  of  loading,  geometry,  and  the  rolling  or  sliding  speed  in  the  line 
or  point  contact.  These  variations  result  in  a  squeeze  effect  which  affects 
the  minimum  film  thickness  distribution.  An  example  of  this  is  the  ball 
bearings  in  a  rotordynami c  system  in  which  there  exist  cyclic  variations  of 
the  dynamic  load.  Recently,  the  transient  hydrodynamic  and  elastohydrodynami c 
line  contact  problem  has  received  much  attention  (Refs.  1  to  3).  Among  the 
several  authors,  Vichard  (Ref.  1)  pioneered  the  basic  transient 
characteristics  of  the  line  contact  problem  analytically  and  experimentally 
including  the  viscous  damping  phenomenon.  In  this  paper,  the  transient 
solution  of  the  hydrodynami cal ly  lubricated  point  contact  presented. 

In  solving  the  point  contact  transient  problem  numerically,  a  fast 
computer  code  is  needed  to  solve  the  two  dimensional  Reynolds  equation  for 
many  time  steps.  Numerical  methods  for  solving  the  simultaneous  equations 
resulting  from  the  d' scretization  of  the  Reynolds  equation  are  usually 
performed  using  either  iterative  methods  or  semidirect  methods  (Ref.  4).  The 
former  commonly  involves  the  Gauss-Seidel  method,  the  latter  combines  the 
Newton-Raphson  method  with  a  direct  inversion  of  the  Jacobian  matrix.  An 
Important  difference  between  the  iterative  method  and  the  semidirect  method  Is 
that  the  Initial  guess  plays  an  Important  role  in  the  latter,  whereas  the 
former  is  relatively  insensitive  to  the  Initial  guess.  With  the  semidirect 
method,  the  use  of  a  previous  solution  as  an  initial  guess  accelerates  the 
solution  process,  but  a  good  initial  guess  usually  does  not  help  the  iterative 
method  significantly  (Ref.  4).  The  semidirect  method  is  preferred  for 
transient  problem  since  the  solution  of  the  previous  time  step  accelerates  the 
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next  step  solution.  Furthermore,  the  Newton-Raphson  method  has  a  quadratic 
convergence  rate,  so,  in  general,  the  solution  can  be  terminated  within  ten 
iterations.  When  a  parallel  processing  computer  using  vectorization  is 
employed  the  matrix  inversion  is  very  fast.  In  addition,  there  is  no  need  to 
use  underrelaxation  factors,  and  the  solution  can  be  obtained  more  rigorously 
than  is  typical  with  iterative  methods.  The  matrix  inversion  can  be  done  by 
the  Thomas  algorithm,  arid  there  is  no  need  to  store  the  whole  Jacobian  matrix. 

When  the  semidirect  method  is  used  in  the  point  contact  problem,  the 
cavitation  boundary,  where  the  Reynolds  boundary  condition  (B.C.)  is  specified, 
is  difficult  to  locate.  There  is  a  fundamental  difference  between  the  line 
contact  and  the  point  contact  problem.  In  the  line  contact  case,  the  Reynolds 
equation  is  integrated  once;  the  Neumann  condition  is  Introduced;  and  the 
Integration  constant  Is  found  as  a  part  of  the  solution.  In  the 
two-dimensional  problem,  the  Reynolds  equation  can  not  be  Integrated.  Since 
the  Reynolds  B.C.  Insures  mass  conservation  across  the  boundary,  the 
cavitation  boundary  should  be  located  as  accurately  as  possible.  However,  the 
location  is  not  known  in  advance;  it  is  a  part  of  the  solution.  It  Is  a  free 
boundary  where  two  B.C.'s  are  present:  Dirichlet  B.C.  (pressure  is  zero),  and 
Neumann  B.C.  (normal  pressure  gradient  is  zero).  The  relaxation  method  of 
Chri stopherson  (Ref.  5),  deriv -u  tor  the  hydrodynamic  lubrication  of  a  journal 
bearing,  has  been  used  to  solve  this  kind  of  free  boundary  value  problem. 

This  method  truncates  negative  computed  pressures  whenever  they  occur  during 
Iteration.  However,  this  method  can  not  be  used  In  the  semidirect  method.  In 
this  work  a  body-fitted  coordinate  system  is  Introduced  which  transforms  the 
unknown  boundary  Into  a  fixed  boundary  and  the  unknown  boundary  function  is 
Introduced  into  the  equations  of  motion.  The  smooth  cavitation  boundary  is 
found  up  to  truncation  and  machine  errors,  whereas  the  result  for 
Chri stopherson' s  method  is  dependent  upon  the  mesh  size  near  the  boundary.  To 
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detect  the  minute  change  of  the  cavitation  boundary  between  the  adjacent  time 
steps,  the  current  method  is  desirable.  Another  advantage  of  this  method  is 
that  a  nonzero  pressure  gradient  condition  can  be  implemented  for  very  lightly 
loaded  cases  where  surface  tension  may  play  an  important  role,  or  for 
non-Newtonian,  viscoelastic  fluids. 

In  the  present  paper  the  transient  hydrodynamic  lubrication  of  a  step 
bearing  is  solved  analytically  to  provide  physical  insight  into  the  transient 
characteristics  of  hydrodynamic  lubrication.  Next,  the  point  contact  problem 
is  solved  numerically  by  the  Newton-Raphson  method  with  Thomas  algorithm. 

This  method  is  fast  and  does  not  require  vast  computer  storage.  Parallel 
processing  by  vectorization  is  also  utilized. 

The  variation  with  time  over  a  loading  cycle  of  the  minimum  film 
thickness,  squeeze  velocity,  and  the  cavitation  boundary  is  studied  for  a  wide 
range  of  excitation  frequencies. 

NOMENCLATURE 

F  dimensionless  load 

Fg  dimensionless  mean  load 

F  right  hand  side  equation  of  discretized  equation 

f  load,  N  (point  contact),  N/m  (step  bearing) 

fg  mean  load,  N  (point  contact),  N/m  (step  bearing) 

G  dimensionless  cavitation  boundary  function 

G'  first  derivative  of  G  with  respect  to  Y 

G"  second  derivative  of  G  with  respect  to  Y 

g  cavitation  boundary  curve  function 

H  dimensionless  film  thickness 

Hg  dimensionless  minimum  film  thickness 

fig  normalized  dimensionless  minimum  film  thickness,  Hg/Hgm 
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dimensionless  minimum  film  thickness  for  Fq 

f 1 1m  thickness,  m 

minimum  film  thickness,  m 

number  of  iteration  of  Newton-Raphson  method 

length  of  the  step  bearing,  m 

reference  length  for  order-of-magni tude  analysis,  m 

number  of  grid  in  £  direction 

number  of  grid  in  n  direction 

normal  direction  vector 

dimensionless  pressure 

pressure,  N/m2 

radius  of  sphere,  m 

residual  vector  of  discretized  equation 
time,  sec 

solution  vector  of  the  discretized  equations 
average  surface  velocity  In  x-direction,  m/sec 
reference  velocity  for  order-of-magni tude-analysi s ,  m/sec 
dimensionless  coordinate  along  rolling  direction 
dimensionless  inlet  boundary  location  in  X-direction 
coordinate  along  rolling  direction 
Inlet  boundary  location  In  x-direction 
dimensionless  coordinate  transverse  to  rolling  direction 
dimensionless  Inlet  boundary  location  in  Y-dlrection 
coordinate  transverse  to  rolling  direction 
inlet  boundary  location  in  y-dlrection 
viscosity-pressure  coefficient,  m2/N 
dimensionless  viscosity-pressure  coefficient 
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fractional  loading  amplitude  for  sinusoidal  loading 
Y  dimensionless  frequency 

5  dimensionless  film  thickness  of  the  step  bearing 

6  normalized  film  thickness  of  step  bearing,  $/Som 

6om  dimensionless  film  thickness  of  the  step  bearing  for  mean  load 

li  lubricant  viscosity.  Pa-sec 

p  dimensionless  lubricant  viscosity 

PO  lubricant  viscosity  at  atmospheric  pressure,  Pa-sec 

v  kinematic  viscosity,  m^/sec 

E,n  coordinates  of  transformed  domain 

p  lubricant  density,  kg/rn^ 

x  dimensionless  time 

4>s  phase  angle  of  the  step  bearing  solution,  deg 

<j>p  phase  angle  of  the  point  contact  solution,  deg 
Q  physical  domain 

Q'  computational  domain 

co  frequency  of  sinusoidal  loading,  ':ycle)/sec 

ANALYTICAL  SOLUTION  OF  A  STEP  BEARING 
Consider  the  simple  step  bearing  shown  in  Fig.  1.  Note  that  the  step 
bearing  used  here  is  subjected  to  an  oscillating  normal  motion  and  is  closed 
at  the  exit  end.  To  the  authors'  knowledge,  this  particular  solution  is  not 
available  in  the  literature  and  is  therefore  presented  here.  The  film  profile 
and  the  dynamic  force  are: 

h< x , t)  =  h(t) ,  0  <  X  <  L, 


—i 

H 

X 

O 

(1) 

f0<l  +  3  sin  cot). 

(2) 
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For  an  Incompressible,  Isovlscous,  Newtonian  fluid,  the  governing 
equation  is, 


h  (h3  I?)  = 


3h  ,,  3h  . 
m  9x  +  2}i0  3x  ' 


U1 

V  T 


The  boundary  conditions  and  the  initial  condition  are, 

\ 

p  =  0  at  x  =  0, 

h  =  0  at  x  *  L, 

h  =  h^  when  t  =  0. 
With  the  following  definitions, 
tu. 


&  -  *■ 
S  -  L- 


m.  y  x . 

t  =  L  ,  a  =  L, 


p  , 

Vm' 


F  = 


F*0U 


m 


the  dimensionless  equations  are. 


8_ 

3X 


(«3  I)  -  '2  '2  Mr. 


ax  T  ,fc  ax' 

S<X,x)  =  S(x) ,  0  <  X  <1 , 

=0,  X  =  1 , 

P  =  0  at  X  =  0, 

6  =  6  j  when  x  =  0, 

F(x)  =  Fgd  +  |3  sin  yx)  • 

After  integrating  Eq.  (5)  three  times  using, 

,1 


J, 


P  dX  =  F(t) 


Y  = 


ojL 


the  following  nonlinear  first  order  differential  equation  results: 

2  96  3  lc,  , 

?  5T-  -  IFU>- 

The  solution  of  Eq.  (8)  subject  to  the  B.C.'s  in  Eq.  (5)  is, 


(3) 


(4) 


(5) 


(6) 
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S(x) 


1 


Fq  fo  &y 

LVS2  6  +  2<y2  +  9) 


.  -3x  0 

le  +  r  + 


(8) 


,-1/2 


2(y  +  9) 


(3  Sin  yx  -  y  COS  yx) 


(9) 


After  a  sufficiently  long  time  <x  ■+  ®),  the  exponential  term  vanishes  to  zero, 
and  the  time  variation  of  the  film  thickness  becomes, 
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Vt 


slnCyt  -  <^  ) 


4»s  =  tan  y. 


The  formula  for  the  squeeze  velocity  Is  obtained  by  differentiating 


Eq.  (10), 


-  (-  8 
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ANALYTICAL  FORMULATION  OF  THE  POINT  CONTACT  PROBLEM 
The  physical  model  Is  Illustrated  In  Fig.  2.  The  radius  of  the  sphere  is 
R  and  the  dynamic  force  is  the  same  as  that  of  the  step  bearing.  The  two 
dimensional,  transient,  Incompressible  form  of  the  Reynolds  equation  for 


Newtonian  flow  is. 


Si  £!*)♦!?(£  8)-'^  8-”  8- 


where 


p  =  p(x,y,t)  I  (12 

h  =  h(x,y,t)  ) 

M  =  p(x,y,t>.  ] 

The  parabolic  approximation  of  the  film  thickness  equation  of  the  sphere 


h  -  hQ  +  y  (x2+  y2).  (1 

At  a  given  time,  the  generated  pressure  distribution  is  balanced  by  the 


dynamic  load, 


P(x,y, t)dx  dy. 
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The  plezovlscous  effect  Is  modelled  by  the  Barus  relation  (16): 


»  =  P0e 


The  boundary  conditions  are: 


p  =  0  at  x  =  x. 


o  <  y  <  yR. 


p  =  0  at  x.  <  x  <  g(yR,t)  y  =  yR( 


p  =  0;  =  0  at  x=  9^y.t)  0  <  y  <  yR,  / 

9n  i 

-  0  at  XA  <  x  <  g(0,t>  y  -  0.  I 

At  the  cavitation  boundary,  x  =  g(y,t),  the  pressure  and  the  normal 
pressure  gradient  are  zero  (Reynolds  B.C.).  Using  symmetry  at  y  =  0,  the 
Neumann  B.C.  Is  Imposed  and  only  half  of  the  domain  Is  modelled. 

With  the  following  definitions, 


X  -  y  _  1-  H  - 

A  —  n  i  1  “  n  >  n  -  n  > 


p  =  -BB_. 

¥m' 


ViriR' 


Y  ■  r; 


the  dimensionless  equations  are: 


3_(H^3P  3_  9P  _  i?  0H  1P  8H 

ax  Vu  a xy  +  ay  Va  ay/  =  u  ax  +  u  dz’ 


H  =  H0  +  ~  (X2  +  Y2), 


rYB  rG(Y) 


F(x)  =  2 


P(X,Y,x)dX  dY, 


0  JX, 


F(i)  =  Fn(  1  +  (3  si n  y  x) , 


M  =  e 


with  a  *  1.5131X10"°  In  this  study. 
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To  fix  the  unknown  cavitation  boundary,  the  following  body-fitted 
coordinate  transformation  shown  in  Fig.  3  is  introduced: 

yb(X  -  XA) 

^  =  G(Y,x>  -  XA’ 

n  -  Y,  (If 

|J|  =  <G<Y.t>  -  Xa)/Yb. 

|J|  is  the  Jacobian  of  the  coordinate  transformation  which  shows  that  as 
long  as  G(Y.-c)  Is  not  equal  to  X/\,  there  exists  a  conformal  mapping  between 
the  physical  domain  and  the  computational  domain. 

The  differentiations  transforms  to  the  following: 

Y 

3_  _  B  8_ 

8X  “  G  -  XA  3£’ 

i_  9_  EG'  _d_ 

8Y  =  3n  "  G  -  XA  9C 


2  Y  7 

<r _ 'b  jr 

9  2  2' 

ax^  <g  -  xAr  ar 


a2  _  s2(G')2  a2  _  2fg*  a2 

3Y2  (G  -  XA)2  as2  G  "  xa  8^9n 


2 ( G ' )  -  G' ' <G-Xa>J  g 
(G  -  XA>2  ^ 


The  Reynolds  equation  in  the  (E,n>  system  is, 


where 


AP__  +  BP-  +  CP  +  DP-  +  EP  +  F  =  0, 
ZZ  En  nn  Z  n 


A  =  A3[y2  +  ^(G1)2], 


B  =  A3C-25G’(G-Xa)]. 


C  =  a3<g  -  xA>fc, 


D  =  A] Yg(G  -  XA>  -  A2CG'(G  -  XA)  +  A3c[2(G')2  -  G"(G  -  XA>]  , 
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E  =  A2<G  -  XA)t, 


F  =  -A4(G  -  XA)  , 


.  3H2  3H  u3  3  M 

Ai  =  —  ax  +  H  ax  -  • 


A  3T  3H  u3  a_  1 

a2  -  3Y  3Y  - 


A  -  — 


,  ,,  3H  3H 

A4  =  12  3X  +  12  3x ' 


In  the  above  formulation,  A] ,  A2,  A3,  and  A4  can  be  transformed  to  the 
(E,n)  coordinate  system  using  Eq.  (19).  At  the  cavitation  boundary, 


1  +  (G1) 


Since  3P/3n  *  0  at  E  =  V 


0,  P  =  0  at  £  =  Y 


At  n  =  0(Y  =  0),  the  symmetry  condition  is, 


3P  3P  EG'  3P  n 
3Y  =  3n  ‘  G  -  XA  3E 


But,  G '  =  0  due  to  the  symmetry  of  cavitation  boundary  and  it  follows 


that, 


=  0  at  n  =  0.  < 

on 

The  transformed  film  thickness  equation  and  the  force  balance  equation 
are  expressed, 


(25) 

(26) 


In  the  above  formulation,  the  unknown  boundary  curve  function  G  is 
introduced  into  the  governing  equations  while  the  computational  domain  is 
f i xed . 

NUMERICAL  METHODS 

Equation  (20)  is  a  nonlinear  partial  differential  equation.  The 
nonlinearity  Is  due  to  the  piezoviscous  relation  and  to  the  function  G  in 
the  transformed  Reynolds  equation. 

Spatial  Discretization 

In  order  to  minimize  the  number  of  grid  points  while  maintaining  accuracy, 
a  smoothly  varying  nonuniform  spacing  is  generated  by  a  two-sided  stretching 
function,  (hyperbolic  tangent)  (Ref.  10).  The  finest  spacing  is  near  the 
cavitation  boundary  which  is  also  near  the  maximum  pressure  gradient. 

Figure  4  shows  the  finite  difference  mesh  structure.  The  increments  in 
£  and  n  and  are  such  that 

*1  '  *1-1  "  A* 

*I+1  -  *1  = 

nj  -  nJ_1  =  An 

nJ+l  -  =  rnAn' 

By  the  Taylor  series  expansion,  the  finite  difference  approximations  of 
derivatives  with  respect  to  £  and  n  are, 
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3P  _  r5PI-l  ,J  *  (rE  "  1)PI,J  +  PI  +  1  ,J 


r^U  +  r^)AE 


on  1  1  +  ~  1>PT  1  +  PT  1  1 

3 P  _  n  I ,  j- 1  n  I , j  I , J+ 1 

8rj  r  ( 1  +  r  )An 

n  n 


02P  o  r£PI-1  ,J  +  1  )PI,J  +  PI  +  1  ,3 


=  2 


3C 


r^(  1  +  r^)AE‘ 


82P  =  2  rnPI,J-1  ~  (rn  *  1)PIJ  *  PI,J-t-l 
3n2  r  (1  +  r  )An2 


2 

3  P 


3^3n  r_r  (1  +  rr)(l  +  r  )A£Ari 
5  n  Z  n 


r2r2P r  .  .  .  -  r2(r2  -  1)PT  .  .  -  r2PT  ,  .  .  -  r2(r2  -  1)PT  .  .  +  (r2  -  1) 
K  n  i-l.J-l  n  Z  I.J-i  n  I+l,j-l  Z  n  I - 1 , J  E 


x  (r2  -  1)PT  +  (r2  -  1)PT  .  ,  -  r2PT  ,  ,  .  +  (r2  -  1)PT  .  .  +  P.  .  .  . 
n  I,J  n  1  +  1  ,J  E  I-l.J+l  E  I, J+l  1+1 ,J+1 


(28) 

Substituting  Eq.  (28)  into  the  transformed  Reynolds  Eq.  (20)  the  following 
discretized  equations  results, 

RI,J  =  C1PI-1,J+1  +  C2PI ,  J+1  +  C3PI  +  1 ,  J+l  +  C4PI-1,J  +  C5PI,J  +  C6PI+1,J 

+  C7  +  C  P  +  C  P  +  C  =  0,  (29) 

7I-1,J-1  8  I.J-l  9  I+l.J-l  10 


wi  th 


P1 ,J  =  PNI,J  =  0  1  <  J  <  NO, 

PI,NJ  =  0  1  <  I  <  NI, 

PI,0  =  PI  ,2  1  <  I  <  NI, 


PNI  +  1  ,J  "  PN I  —  1  ,J  1  <  J  <  NJ. 
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Steady-State  Solution  Method 

The  transient  solution  is  formed  by  computing  the  steady-state  solution 
for  each  time  step  including  the  squeeze  term.  The  numerical  technique  for 
the  steady-state  solution  along  with  the  Thomas  algorithm  and  Newton-Raphson 
method  is  described  first. 

The  discretized  form  of  transformed  equation  is, 

K(u)u  =  F  (30) 

The  vector  u  represents  the  unknown  values,  pressures  and  cavitation 
boundary.  For  an  isoviscous  condition  K(u),  contains  the  function  G,  and, 
for  a  piezoviscous  condition,  it  includes  pressures  as  well.  The  discretized 
simultaneous  equations  are  nonlinear.  Even  for  the  linear  free  boundary  value 
problem,  it  has  a  nonlinear  characteristics  since  the  unknown  boundary  is 
associated  with  the  solution. 

The  Newton-Raphson  method  is  described, 

Vi  -  \  -  r,<V«<V  (31> 

where  R(uk)  =  K< u^) u^  -  F  is  the  residual  vector  and  J(u)  is  the  Jacobian 
of  the  system  of  equations.  In  practice,  the  iteration  is  organized  as, 

J(uk)Auk  =  -R(uk>,  uk+1  =  uk  +  Auk.  (32) 


For  this  study,  the  vector  u  is, 

u  "  (P2 , J ’  P3,J . PN I — 1 , J ,Gj)T’  J  =  1  ,NJ  -  1  (33) 

in  which  Pi  j  and  Pni.j  are  zero  from  the  Dlrichlet  boundary  condition. 

The  residual  vector  R  is, 


R 


rni-i  ,j 


J  =  1 ,  NJ  -  1 


(34) 
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The  Jacobian  matrix  is  a  block  tridiagonal  matrix  in  Fig.  5,  and  each 
block  is  a  one-sided  arrow-shaped  matrix,  Fig.  6.  In  the  formulation  of  each 
block  matrix  of  the  Jacobian,  the  last  columns  are  the  di fferentiations  of  the 
residual  vector  with  respect  to  the  cavitation  boundary  function,  G.  Since  all 
the  coefficients  in  the  discretized  Reynolds  equation  are  composed  of  Gj,  G'j, 
and  G''j,  it  is  easier  to  calculate  them  numerically  (Ref.  11)  using: 


3R 


I,J  1 


9G- 


4,j(Gj  +  cg,wi,j)  “  ri,j(Gj,wi,j) 


(35) 


where  WI,J  contains  all  other  variables  except  Gj.  The  value  of 

eg  can  be  chosen  to  be  sufficiently  small  not  only  to  maintain  good  accuracy 

of  Eq.  (35)  but  to  prevent  serious  round-off  errors.  In  this  calculation, 
eg  is  set  to  10-9  in  double  precision. 

The  block  tridiagonal  system  of  Eq.  (30)  is  solved  by  the  Thomas 

algorithm  (Ref.  12).  This  algorithm  inverts  the  whole  matrix  at  a  time  by 

matrix  multiplication  and  inversion  of  the  block  matrix,  which  is  quite  fast 
on  a  parallel  processing  computer  with  small  memory  storage  size  equal  to 
2  x  NI  x  NJ  x  NJ.  The  matrix  inversion  is  accomplished  using  UNPACK. 

The  Newton-Raphson  method  requires  a  good  initial  guess  of  the  solution. 
For  this  purpose,  the  Gauss-Seidel  iteration  method  is  used  to  get  an 
approximate  pressure  distribution  and  cavitation  boundary  location.  Once  one 
solution  is  obtained  by  the  Newton-Raphson  method,  it  is  used  for  the  guess  to 
next  solution.  The  convergence  criteria  are 


(1)  pressure 


1 .0x1  O'4 
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(2)  cavitation  boundary 


1 .0x10  4 


(3)  force  balance 

F  input  ~  ^output  . 

- P -  <  1.0x10"^ 

r i nput 

In  order  to  make  sure  of  the  convergence,  the  L.2-norm  of  the  residual 
vector  is  monitored.  In  general,  the  solution  converges  within  3  to  8 
iterations.  In  this  study,  NI  =41,  NJ  =  31 . 

Transient  Solution  Method 


For  the  steady-state  solution,  the  problem  is  to  find  Hq  for  a  given 
load,  or  for  a  hydrodynamic  case,  the  load  capacity  can  be  calculated  for  a 
given  Hq.  But,  for  the  transient  case,  there  is  an  additional  unknown  value 
to  be  determined,  the  squeeze  velocity.  The  basic  solution  technique  is  to 
use  a  "time-march."  That  is,  Hq  is  fixed  from  the  previous  time  step,  and 
the  squeeze  velocity  is  found  that  balances  the  generated  pressure  distribution 
with  the  dynamic  force  at  that  time.  The  detailed  computation  procedure  is 
provided  in  Fig.  7.  At  the  first  time  step,  the  steady-state  Reynolds 
equation  is  solved  to  find  HQm-  and,  fixing  Hq,  the  transient  Reynolds 
equation  is  solved  including  the  squeeze  term  to  find  the  squeeze  velocity 
using  the  force  balance  equation.  For  this  purpose,  a  bisection  method  is 
used,  with  an  approximate  range  of  squeeze  velocities  according  to  the  history 
of  dynamic  force  and  the  minimum  film  thickness  variation.  Once  a  converged 
solution  is  obtained,  the  minimum  film  thickness  of  the  next  time  step  is 
estimated  from  the  following  expression: 
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(36) 


n 

At, 


(n  =  present  time  step) 


The  film  thickness  and  squeeze  velocities  are  established  at  successive 
time  steps  and  the  calculation  is  continued  beyond  the  first  complete  loading 
cycle  until  the  periodic  requirement  is  reached.  The  convergence  criterion  is: 


(H0n  "  Vn 


IC 

(H0n 


_4 

<  1.0x10  ,  IC  =  number  of  cycle 


In  this  calculation,  361  time  steps  with  1°  increment  are  used  in  one 
loading  cycle. 


RESULTS  AND  DISCUSSION 


The  analytical  solution  of  the  step  bearing  demonstrates  that  <5  approaches 
one  with  a  phase-lag  of  90°  as  y  increases  (Fig.  8).  This  asymptotic 
behavior  is  due  to  the  squeeze  action  caused  by  the  dynamic  forces.  Figure  9 
shows  the  squeeze  variation  of  Eq.  (11).  This  phenomenon  is  physically 
similar  to  a  nonlinear  massless  spring-damper  system  with  forced  vibration 
shown  in  Fig.  10,  sometimes  referred  as  a  "half  a  degree  of  freedom  system." 

The  response  of  this  system  is  that  the  amplitude  approaches  a  constant  value 
and  the  phase-lag  goes  to  90°.  Although  the  transient  solution  of  the  point 
contact  problem  can  not  be  solved  analytically  and  requires  numerical 
computation,  it  may  be  speculated  that  basically  it  also  has  a  similar 
nonlinear  spring-damper  system.  In  the  following  example,  the  numerical 
results  of  the  point  contact  problem  are  compared  to  the  step  bearing  solution. 

For  this  study,  Fq  =  3000  and  3  =  0.3  with  different  y's.  The 
minimum  film  thickness  for  Fq  Is  1.2471x10“ 5  for  the  isoviscous  case  and 
1.3907xl0_5  for  the  piezoviscous  effect  with  =  0.08,  Yb  =  0.06.  Figure  11 
shows  the  pressure  distribution  for  Fq  and  Fig.  12  delineates  the  detailed 
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cavitation  boundary  curve  in  which  the  minimum  value  of  G  occurs  at  Y  =  0 
and  it  increases  up  to  a  certain  location  and  then  decreases  because  of  the 
geometry  of  the  sphere. 

Figure  13  illustrates  the  time  variation  of  the  normalized  minimum  film 
thickness  <R())  during  one  loading  cycle  with  361  time  steps.  The  squeeze 
velocity  distribution  is  shown  in  Fig.  14  for  different  y's.  These  results 
are  qualitatively  similar  to  those  of  the  step  bearing  solution.  However,  it 
should  be  noted  that  the  order  of  the  nondimens ional  excitation  frequencies  is 
different  since  L  is  used  as  a  reference  length  in  the  step  bearing  while 
R  is  used  for  the  point  contact  case. 

Equation  (10)  may  be  put  in  the  following  form. 


S  = 


1  +  as3  sin(yx  - 


1/2 


where 


as  = 


4>s  =  tan'  (Xsy> 


xs  =  3* 


1  ♦  (Xsy>‘ 


(37) 


The  variation  of  as  and  4>s  are  plotted  in  Fig.  15. 

For  quantitative  analysis  of  the  transient  point  contact  problem,  the 
following  formula  is  suggested  by  Eq.  (37), 


ap|3  s i n ( Yt  -  $p>. 


(38) 


Equation  (38)  is  deduced  based  on  the  fact  that  Hom  is  Inversely 
proportional  to  Fq  whereas  SQm  to  ^Fq .  The  unknown  values  in  Eq.  (38), 

ap  and  4>p,  are  obtained  by  a  nonlinear  least  square  fit  with  361  data 
point.  Figure  16  shows  the  comparison  between  the  numerical  results  and  the 
curve  fit.  The  best  curve  fit  can  be  obtained  by  letting  the  numerator  of 
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Eq.  (38)  be  variable,  however,  it  is  near  1,  for  example,  1.005  for  y  =  100, 
1.019  for  y  =  1000.  The  curve  fitting  results  are  recorded  in  Table  1. 

Fig.  17  shows  the  variation  of  ap  and  <j>p,  qualitatively,  similar 
characteristics  to  the  analytical  step  bearing  solution  with  different  order 
of  magnitude  of  y  (Fig. 15).  The  value  of  Xp  is  obtained  assuming  the 
following  relation, 

4>  =  tan_1(\  v)  (39) 

P  P 

Xp  is  nearly  constant  over  a  wide  range  of  y,  approximately  0.0054.  If  an 
analytical  solution  were  possible,  the  ap  would  be  a  function  of  Xp. 

However,  since  it  also  would  be  a  function  of  the  geometry  associated  with  the 
cavitation  boundary,  no  attempt  is  made  to  obtain  a  form  similar  to  Eq.  (37). 
Instead,  for  design  purposes,  Eq.  (38)  can  be  used  along  with  Table  1. 

For  the  piezovlscous  solution,  ap  is  smaller  than  that  of  the  Isoviscous 
solution  (Fig.  17),  but  4>p's  are  virtually  the  same.  The  ap's 
asymptotically  approach  those  of  the  isoviscous  case.  Figure  18  shows  this 
more  vividly.  Due  to  the  piezoviscous  effect,  the  distribution  of  flo  1s 
more  damped  with  the  same  phase  angle.  The  Xp's  for  the  piezovlscous  case 
are  nearly  constant  and  equal  to  the  isoviscous  case  (see  Table  1).  This 
Implies  that  Xp  is  a  characteristic  of  the  transient  point  contact  problem 
of  the  current  model. 

Figure  19  Illustrates  the  location  of  the  outlet  boundary  at  Y  =  0 
normalized  by  that  for  the  steady-state  solution  of  Fq.  For  the  steady  state 
case,  G(0,t)  approaches  the  point  of  contact  as  the  load  Increases.  However, 
when  y  is  greater  than  zero  there  exists  a  substantial  variation  in  G(0,t) 
due  to  the  squeeze  action.  When  the  squeeze  is  downward,  G(0,x)  may  be 
stretched  outward  and  vice  versa.  For  example,  when  y  =  200,  there  is  a 


19 


downward  action  between  a-b  and  c-d  in  Fig.  19,  and  upward  action  between 
b  and  c.  These  points  correspond  to  those  in  Fig.  14. 

In  the  foregoing  analysis,  the  Reynolds  Eq.  (12)  neglects  the  inertia 
forces.  But,  as  y  Increases,  the  validity  of  this  assumption  becomes 
suspect.  This  assumption  is  examined  by  an  order-of-magni tude-analysi s 
of  the  steady-state  Navier-Stokes  equation  in  Ref.  (12).  When  the  modified 
Reynolds  number  is  much  less  than  one, 


(40) 


the  inertia  forces  can  be  neglected.  Here,  uq  Is  a  reference  velocity, 

Q.  is  a  reference  length  In  the  x-dlrection,  and  hg  is  that  in  the  film 
thickness  direction.  Using, 

UQ  =  Ru;  Y  *  (41) 

m 

the  following  relation  for  the  validity  of  the  assumption  that  inertia 
forces  are  negl igible  is, 


Y 


H0Re- 


Re 


(42) 


For  example,  if  Hq  =  lO-^  r  =  io-2  m,  u  =  0.1  m/sec,  and 
v  =  10" ^  m2/sec, 

y  «  108  (43) 

Even  for  y  =  1000,  inertia  effects  remain  negligible. 


CONCLUSIONS 

The  transient  solution  of  the  hydrodynamical 1 y- 1 ubri cated  point  contact 
problem  including  the  squeeze  effect  is  obtained  numerically  using  the 
ball-on-plane  model.  A  new  computational  algorithm  is  Implemented  to  deal 
with  the  cavitation  boundary  by  the  semldirect  method  with  the  advantage  of 
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supercomputi ng .  This  method  provides  a  faster  and  more  rigorous  way  to  solve 
the  nonconformal  contact  problem  with  a  Newtonian  fluid  than  the  conventional 
iterative  method,  and  the  flexibility  to  deal  with  more  complex  boundary 
conditions  for  lightly  loaded  bearings  and  more  realistic  rheological  models. 

The  qualitative  and  quantitative  analysis  is  compared  with  the  analytical 
solution  of  a  dynamically  loaded  step  bearing  solution  using  a  nonlinear  curve 
fitting  method.  It  is  found  that  there  exists  a  character i sti c  similarity  in 
the  transient  responses  to  a  nonlinear  massless  (i.e.,  no  inertia) 
spring-damper  system,  in  terms  of  the  variation  of  the  minimum  film  thickness 
and  phase  angle.  According  to  an  order  of  magnitude  analysis,  it  is  confirmed 
that  the  inertia-forces  are  negligible  for  a  wide  range  of  practical 
excitation  frequencies. 

These  results  can  be  applied  to  the  design  of  moderately  loaded  ball 
bearings  in  rotordynamic  systems  and  can  be  extended  to  gear  design  adding  the 
time  variation  of  the  geometry  and  speed.  For  highy  loaded  elliptical  contact 
case,  the  elastic  deformations  and  ellipticity  parameter  need  to  be  considered 
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TABLE  1  -  CURVE  FITTING  RESULTS  OF  EQ.  (38) 


Isovi scous 

Piezovi scous 

Y 

aP 

XP 

mm 

o 

XP 

0 

0.890 

0.0 

0.734 

0.0 

25 

.882 

7.0 

0.00491 
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0.00491 
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.00523 
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.00526 
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28.3 
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28.4 
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.00544 
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47.4 

.00544 
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53.7 
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.309 
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.00541 
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69.8 
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69.8 

.00544 
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76.1 

.00539 

.176 

76.1 

.00539 

1000 

.159 

79.4 

.00536 

.134 

79.4 

.00536 
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Figure  1 .  -  Schematic  view  ol  the  step  bearing  configuration. 
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F  gure  6  •  The  elements  in  diagonal  Wock  matrix 


Figure  7.  •  Flow  digram  lor  the  transient  solution. 
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STEADY  STATE. 


TIME  STEP 

TIME  STEP 

Figure  8.  -  Normalized  film  thickness  versus  nondimensional 
time  during  one  loading  cyde  for  the  step  bearing. 

Figure  9  ■  Squeeze  velocity  divided  by  versus 

nondimensional  bme  during  one  loading  cycle 
for  the  step  bearing. 

F0  0.3000E  04 

MAXIMUM  P:  0  2809E  08  A 

HO  0.1247E-04  Al 

■  T  ^VXVNN 


Figure  10.  -  HaJf  a  degree  of  freedom  model. 


Figure  11-  The  pressure  distribution  of  half  of  the  domain  for  the 
point  contact  problem  for  F0  *  3000 
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16  Abstract 

The  transient  analysis  of  hydrodynamic  lubrication  of  a  point-contact  is  presented.  A  body-fitted  coordinate 
system  is  introduced  to  transform  the  physical  domain  to  a  rectangular  computational  domain,  enabling  the  use  of 
the  Newton-Raphson  method  for  determining  pressures  and  locating  the  cavitation  boundary,  where  the  Reynolds 
boundary  condition  is  specified.  In  order  to  obtain  the  transient  solution,  an  explicit  Euler  method  is  used  to 
effect  a  time  march.  The  transient  dynamic  load  is  a  sinusoidal  function  of  time  with  frequency,  fractional 
loading,  and  mean  load  as  parameters.  Results  include  the  variation  of  the  minimum  film  thickness  and  phase-lag 
with  time  as  functions  of  excitation  frequency.  The  results  are  compared  with  the  analytic  solution  to  the  transient 
step  bearing  problem  with  the  same  dynamic  loading  function.  The  similarities  of  the  results  suggest  an 
approximate  model  of  the  point  contact  minimum  film  thickness  solution.. 
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